###############################
# Analysis of conjoint data 
###############################

rm(list=ls())

# load packages 
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)

sink("~/Dropbox/Atacama 2015/10_replication/002_conjoint_analysis.txt")

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

##############################
# Prepare data
##############################

# load data
d = read.dta("~/Dropbox/Atacama 2015/10_replication/conjoint_paipote_final.dta")
d_match = read.dta("~/Dropbox/Atacama 2015/10_replication/conjoint_paipote_matched_sample_final.dta")

dim(d)
dim(d_match)

# drop failed conjoints
d = d[d$failcon!=1,]

# gen treatment indicator
d$t_ind = NA
d$t_ind[d$a2>2]=1
d$t_ind[d$a2<3]=0
table(d$t_ind)

# gen control indicator
d$c_ind = 1 - d$t_ind 
table(d$c_ind)

# gen affected area indicator 
table(d$zone)
d$zone2 = NA
d$zone2[d$zone==1] = 1
d$zone2[d$zone==3] = 1
d$zone2[d$zone==2] = 0
table(d$zone2)

######################################################
## ACME using regression, center as reference category
######################################################

# Check levels
print(levels(d$atideology))  
d$atideology = factor(d$atideology,levels(d$atideology)[c(2,1,4,3)])
print(levels(d$atideology)) 

print(levels(d_match$atideology))  
d_match$atideology = factor(d_match$atideology,levels(d_match$atideology)[c(2,1,4,3)])
print(levels(d_match$atideology)) 

# Heterogenous AMCE: all attributes, treatment indicator, and interactions 
lmout1_center <- lm(selected ~ atideology + atprofession + atgender + atage + atexperience + atexpectations + 
                        t_ind + 
                        atideology*t_ind + atprofession*t_ind + atgender*t_ind + atage*t_ind + atexperience*t_ind + atexpectations*t_ind,
                        data=d)
lmout1_center_c <- round(coeftest(lmout1_center, vcov = vcovCluster(lmout1_center, cluster = d$idnum)),4)
lmout1_center_c

# Heterogenous AMCE: all attributes, control indicator and interactions
lmout2_center <- lm(selected ~ atideology + atprofession + atgender + atage + atexperience + atexpectations + 
                        c_ind + 
                        atideology*c_ind + atprofession*c_ind + atgender*c_ind + atage*c_ind + atexperience*c_ind + atexpectations*c_ind,
                        data=d)
lmout2_center_c <- round(coeftest(lmout2_center, vcov = vcovCluster(lmout2_center, cluster = d$idnum)),4)
lmout2_center_c

###########################
# Multiple Comparisons
###########################

# BH correction
p_atideologyRight_t_ind = lmout1_center_c[14,4]
p_atideologyLeft_t_ind = lmout1_center_c[15,4]
p_atideologyIndependent_t_ind = lmout1_center_c[16,4]

p_no_correction =  c(p_atideologyRight_t_ind,
                     p_atideologyLeft_t_ind,
                     p_atideologyIndependent_t_ind)
  
p_correction_bh =  p.adjust(p_no_correction, "BH")
p_correction_bon =  p.adjust(p_no_correction, "bonferroni")

colnames = c("Right*Treatment",
             "Left*Treatment",
             "Independent*Treatment")

p_no_correction  = round(p_no_correction,4)
p_correction_bh  = round(p_correction_bh,4)
p_correction_bon = round(p_correction_bon,4)

table = data.frame(colnames,p_no_correction,p_correction_bh,p_correction_bon)
colnames(table) <- c("Attributes", "No correction", "BH correction","Bon. correction")
table

stargazer(table, 
          type = "text",
          colnames = TRUE,
          summary = FALSE,
          title="P-values before and after multiple comparison correction", 
          digits=4, 
          rownames=FALSE, 
          float = TRUE, 
          float.env = "table", 
          table.placement = "H")

#############################################
# Robustness check
#############################################

# Different treatment indicator 
rob_lmout1_center <- lm(selected ~ atideology + atprofession + atgender + atage + atexperience + atexpectations + 
                                  zone2 + 
                                  atideology*zone2 + atprofession*zone2 + atgender*zone2 + atage*zone2 + 
                                  atexperience*zone2 + atexpectations*zone2,
                                  data=d)
rob_lmout1_center_c <- round(coeftest(rob_lmout1_center, vcov = vcovCluster(rob_lmout1_center, cluster = d$idnum)),4)
rob_lmout1_center_c

# Matched sample
d_match$t_ind = NA
d_match$t_ind[d_match$a2>2]=1
d_match$t_ind[d_match$a2<3]=0
table(d_match$t_ind)

rob_lmout2_center <- lm(selected ~ atideology + atprofession + atgender + atage + atexperience + atexpectations + 
                                  t_ind + 
                                  atideology*t_ind + atprofession*t_ind + atgender*t_ind + atage*t_ind + atexperience*t_ind + atexpectations*t_ind,
                                  data=d_match)
rob_lmout2_center_c <- round(coeftest(rob_lmout2_center, vcov = vcovCluster(rob_lmout2_center, cluster = d_match$idnum)),4)
rob_lmout2_center_c

#####################
# Prepare plots
#####################

# Unexposed respondents and center as baseline
lmout1_center_c 
coeff = lmout1_center_c[,1] 
coeff2 = cbind(coeff) 
coeff3 = coeff2[1:12,] 
c_se = lmout1_center_c[,2] 
c_se2 = cbind(c_se) 
c_se3 = c_se2[1:12,] 
results = data.frame(coeff3,c_se3) 
colnames(results) = c("y1","r1") 
results = results[-1,] 
baseline =  c(0,0) 
results <- rbind(baseline,results[1:3, ], baseline,results[4:5, ],baseline,results[6, ],
                 baseline,results[7:8, ],baseline,results[9:10, ],baseline,results[11, ]) 
rownames(results) = c("1b.atideology","2.atideology","3.atideology","4.atideology",
                      "1b.atprofession","2.atprofession","3.atprofession",
                      "1b.atgender","2.atgender", 
                      "1b.atage","2.atage","3.atage",
                      "1b.atexperience","2.atexperience","3.atexperience",
                      "1b.atexpectations","2.atexpectations")
results
write.table(results, "~/Dropbox/Atacama 2015/10_replication/center_nonexposedgroup.txt", sep="\t")

# Exposed respondents and center as baseline
lmout2_center_c
coeff = lmout2_center_c[,1] 
coeff2 = cbind(coeff) 
coeff3 = coeff2[1:12,] 
c_se = lmout2_center_c[,2] 
c_se2 = cbind(c_se) 
c_se3 = c_se2[1:12,] 
results = data.frame(coeff3,c_se3) 
colnames(results) = c("y1","r1") 
results = results[-1,] 
baseline =  c(0,0)
results <- rbind(baseline,results[1:3, ], baseline,results[4:5, ],baseline,results[6, ],
                 baseline,results[7:8, ],baseline,results[9:10, ],baseline,results[11, ]) 
rownames(results) = c("1b.atideology","2.atideology","3.atideology","4.atideology",
                      "1b.atprofession","2.atprofession","3.atprofession",
                      "1b.atgender","2.atgender", 
                      "1b.atage","2.atage","3.atage",
                      "1b.atexperience","2.atexperience","3.atexperience",
                      "1b.atexpectations","2.atexpectations")
results
write.table(results, "~/Dropbox/Atacama 2015/10_replication/center_exposedgroup.txt", sep="\t")

# Difference Plot from unexposed to exposed: center as baseline 
lmout1_center_c
coeff = lmout1_center_c[,1] 
coeff2 = cbind(coeff) 
coeff3 = coeff2[14:24,] 
c_se = lmout1_center_c[,2] 
c_se2 = cbind(c_se) 
c_se3 = c_se2[14:24,] 
results = data.frame(coeff3,c_se3) 
colnames(results) = c("y1","r1") 
baseline =  c(0,0) 
results <- rbind(baseline,results[1:3, ], baseline,results[4:5, ],baseline,results[6, ],
                 baseline,results[7:8, ],baseline,results[9:10, ],baseline,results[11, ]) 
rownames(results) = c("1b.atideology","2.atideology","3.atideology","4.atideology",
                      "1b.atprofession","2.atprofession","3.atprofession",
                      "1b.atgender","2.atgender", 
                      "1b.atage","2.atage","3.atage",
                      "1b.atexperience","2.atexperience","3.atexperience",
                      "1b.atexpectations","2.atexpectations")
results
write.table(results, "~/Dropbox/Atacama 2015/10_replication/center_differenceplot.txt", sep="\t")

sink()
